function varargout=SteadyState_Solver(PP,Options)
% clear;
% clc;
% PP          =   Setup_PP(struct());

arguments
    PP struct;
    Options.PSI_0 double            =   1;
    Options.FracHN_w_0 double       =   1;
    Options.FracHN_wN_0 double      =   PP.ExoState.Idio_RI.InvDist(2)/PP.ExoState.Idio_RI.InvDist(1);
    Options.DiscFactor_0 double     =   0.02;
    Options.EqualWageFlag logical   =   1;
    
end

%% Preliminaries
PSI_0       =   Options.PSI_0;
DiscFactor_0=   Options.DiscFactor_0;
FracHN_wN_0 =   Options.FracHN_wN_0;
if ~Options.EqualWageFlag
    FracHN_w_0  =   Options.FracHN_w_0;
end

%% Initiate the SS
if PP.InfoPrint>=1
    if Options.EqualWageFlag
        fprintf(['Start Solving the Steady State w/ Equal Wages ...\n']);
    else
        fprintf(['Start Solving the Steady State w/o Equal Wages ...\n']);
    end
    
end
TimeStart   =   tic;
SS          =   struct();
%--------------------------------------------------------------------------
% Exogenously Fixed
%   Note: 
%--------------------------------------------------------------------------
SS.Pir      =   PP.Pir_SS;
SS.Pir_N    =   SS.Pir;
SS.Pir_H    =   SS.Pir;
SS.Pir_F    =   SS.Pir;
SS.Pir_H_s  =   SS.Pir_H;
SS.Pir_F_s  =   SS.Pir_F;

SS.p_N      =   1;
SS.p_T      =   1;
SS.p_H      =   1;
SS.p_F      =   1;

SS.dE       =   SS.Pir_H-SS.Pir_H_s;

SS.N        =   PP.N_SS;              
SS.ir_dom   =   PP.i_dom_SS;
SS.ir_ext   =   PP.i_ext_SS;

[DomSh,EffReturn,KAPPA,FinReturn] ...
            =   PortfolioAdjCost('r_dom',SS.ir_dom,'r_ext',SS.ir_ext,...
                                 'Share_Target',PP.Fin_DomShare,...
                                 'AdjCost',PP.Fin_AdjCost).OptChoice;
PortfolioMat=   [DomSh.Pos,DomSh.Neg; EffReturn.Pos,EffReturn.Neg;...
                 KAPPA.Pos,KAPPA.Neg; FinReturn.Pos,FinReturn.Neg];
SS.PortfolioInfo=PortfolioMat(:);

SS.ir       =   SS.ir_dom;
SS.ir_s     =   SS.ir_ext;
SS.rr_dom   =   SS.ir-SS.Pir;
SS.rr_ext   =   SS.ir_s+SS.dE-SS.Pir;

SS.TAU      =   PP.TAU_SS;
SS.T        =   PP.T_SS;

[SS.O_StateID,SS.O_TargetTrProbMat] ...
            =   MarkovChain_IndependentMixture({PP.ExoState.Idio_FI.IndNode,PP.ExoState.Idio_FI.TrMat},...
                                               {PP.ExoState.Idio_RI.IndNode,PP.ExoState.Idio_RI.TrMat});
SS.FracT_pC =   PP.Cons_FracT*(SS.p_T)^(1-PP.Cons_ElasTN);
SS.FracN_pC =   (1-PP.Cons_FracT)*(SS.p_N)^(1-PP.Cons_ElasTN);
SS.FracH_pC =   PP.Cons_FracH*(SS.p_H/SS.p_T)^(1-PP.Cons_ElasHF)*SS.FracT_pC;
SS.FracF_pC =   (1-PP.Cons_FracH)*(SS.p_F/SS.p_T)^(1-PP.Cons_ElasHF)*SS.FracT_pC;
%--------------------------------------------------------------------------
% Approximation Structure for Value Function, Policy Function, and Distribution
%   Note: 
%--------------------------------------------------------------------------
[SS.FunApp,SS.PolApp] ...
            =   SteadyState_FunPolApp(PP,SS);
        
SS.O_StateIdx=  zeros(PP.ExoState.IdioInc.N,size(SS.O_StateID,1));
for ii=1:PP.ExoState.IdioInc.N
    for jj=1:size(SS.O_StateID,1)
        SS.O_StateIdx(ii,jj)    =   find(ismember(SS.FunApp.ExoState.State,[ii,SS.O_StateID(jj,:)],'rows'));
    end
end
SS.DistApp  =   SteadyState_DistApp(PP,SS);


%% Solve the SS by Finding to Equilibrium 
%--------------------------------------------------------------------------
% Solve the root by Newton Update 
%   Note: 
%       1. Save time by recycling the UCoef from the intermediate step.
%       2. Fastest, but might not be robust when Diff(PSI,ra) is too flat.
%--------------------------------------------------------------------------
SS_Res_Setup=   struct('Bellman_ValItNum',0,'Bellman_ColItNum',1,...
                     'Bellman_TolItNum',50,'Bellman_DampenCoef',1,...
                     'PrintFlag',PP.InfoPrint,'EqualWageFlag',Options.EqualWageFlag);
                     
if Options.EqualWageFlag
    Input_0     =   [PSI_0;DiscFactor_0;FracHN_wN_0];
else
    Input_0     =   [PSI_0;FracHN_w_0;DiscFactor_0;FracHN_wN_0];
end

[SS,ExitFlag]=  SteadyState_Solver_Newton(PP,SS,SS_Res_Setup,Input_0);
%--------------------------------------------------------------------------
% Recollect the Steady State
%   Note: 
%--------------------------------------------------------------------------
TimeElapsed     =   toc(TimeStart);
if PP.InfoPrint>=1
    fprintf(['Steady State is Solved, Elapse Time=',num2str(TimeElapsed,'%.2f'),'s','\n',...
             'Equilibrium Residual: \n',...
             'Labor Market, Total Labor Quantity: ',num2str(SS.RES(2),'%.2e'),', \n', ...
             'Labor Market, H-N Ratio: ', num2str(SS.RES(3),'%.2e'),', \n', ...
             'Domestic Bond Market: ',num2str(SS.RES(1),'%.2e'),', \n',...
             'Final Good Market: ', num2str(SS.RES_FinalGood),'.\n']);
end

%% State Space Reduction
%--------------------------------------------------------------------------
% Distribution
%--------------------------------------------------------------------------
[SS.Dist,SS.CopulaF] ...
            =   DimReduction_Distribution(PP,SS,'Zip',SS.Dist_fzo.QW_fzo);
%--------------------------------------------------------------------------
% Value Function
%--------------------------------------------------------------------------
SS.ValFun   =   SS.UCoef-SS.UCoef;

%% Important Statistics
SS.AggStat  =   [SS.Agg.Avg.c.dom;SS.Agg.Avg.c.ext;SS.Agg.Avg.c.N;SS.Agg.Avg.c.H;SS.Agg.Avg.c.poor;SS.Agg.Avg.c.rich;...
                 SS.Agg.Avg.c.dom_N;SS.Agg.Avg.c.dom_H;SS.Agg.Avg.c.ext_N;SS.Agg.Avg.c.ext_H;...
                 SS.Agg.Avg.c.dom_poor;SS.Agg.Avg.c.dom_rich;SS.Agg.Avg.c.ext_poor;SS.Agg.Avg.c.ext_rich;...
                 SS.Agg.Avg.c.H_poor;SS.Agg.Avg.c.H_rich;SS.Agg.Avg.c.N_poor;SS.Agg.Avg.c.N_rich];

%% Residuals
UnitNum     =   [SS.DistApp.fzo.UnitNum(1),...
                 PP.ExoState.IdioInc.N,PP.ExoState.Idio_FI.N,PP.ExoState.Idio_RI.N]';
SS.O_AvgTrProbMat ...
            =   TrProbMat2AvgTrProbMat(SS.TrProb.UnitTrMat.ReEval,SS.Dist_fzo.QW_fzo,  ...
                                       UnitNum,[3,4]);
SS.O_AvgTrProbMat_FI ...
            =   TrProbMat2AvgTrProbMat(SS.TrProb.UnitTrMat.ReEval,SS.Dist_fzo.QW_fzo,  ...
                                       UnitNum,[3]);
SS.O_AvgTrProbMat_RI ...
            =   TrProbMat2AvgTrProbMat(SS.TrProb.UnitTrMat.ReEval,SS.Dist_fzo.QW_fzo,  ...
                                       UnitNum,[4]);

if Options.EqualWageFlag
    PP.Cons_FracT   =   SS.Cons_FracT;
    varargout{1}    =   SS;
    varargout{2}    =   PP;
else
    varargout{1}    =   SS;
    varargout{2}    =   PP;
end
